Numerical algorithms of the paper: 
Title:   Fast, Accurate and Robust Adaptive Finite Difference Methods for Fractional Diffusion Equations
Journal:  Numerical Algorithms
Authors:   Santos B. Yuste* and J. Quintana-Murillo  
Departamento de Física and Instituto de Computación Científica Avanzada (ICCAEX), Universidad de Extremadura, 06006 Badajoz, Spain
* santos@unex.es

Mathematica  version...

StringJoin[{"The running Mathematica version is ",ToString[$Version]}]

Numerical algorithms for the T&E and predictive adaptive methods.

We start implementing the  L1 finite difference scheme with variable timesteps for solving the fractional PDE 
	
The finite difference scheme is given by Eq. (10):

 where 

 or, in matrix-vector form, by Eq. (11):

The coefficients     of the paper are given here by the function TT[m,n,γ]  where  (Subscript[t, n])=tt[n]

TT[m_,n_,gamma_]:=((tt[n ]-tt[m])^(1-gamma)-(tt[n ]-tt[m+1])^(1-gamma))/(tt[m+1]-tt[m])

The coefficient
   
of the paper is given here by the function GG[j,n,γ] , where  (Subscript[F, j]^((n)))=Γ(2-γ) ((((Subscript[t, n])-(Subscript[t, n-1])))^(γ ))f((Subscript[x, j,])(Subscript[t, n])) [see Eq. (9)], and f((Subscript[x, j,])(Subscript[t, n]))  is given by  funSource[ j Δx,(Subscript[t, n]),γ] .  The solution U is stored in the matrix Umatrix. The element (1+m,1+j) of this matrix, i.e.,  Umatrix[[1+m,1+j]],  represents  (Subscript[U, j]^((m))) with j=0,1,..., m=0,1,...

GG[j_,1,gamma_]:=Umatrix[[1,j+1]] +
(tt[1]-tt[0])^gamma*Gamma[2-gamma]*funSource[j*deltax,tt[1],gamma]    (* El elemento 1 de Umatrix es la solución en el instante inicial: u(x,0)   *)

GG[j_,n_,gamma_]:=(
kk=(tt[n]-tt[n-1])^gamma*Sum[(Umatrix[[m+2,j+1]]-Umatrix[[m+1,j+1]]) *
TT[m,n,gamma],{m,0,n-2}];
Umatrix[[n,j+1]] -kk+(tt[n]-tt[n-1])^gamma*Gamma[2-gamma]*funSource[j*deltax,tt[n],gamma])

The coefficient (Subscript[S, n]) [Eq. (7)]  is given by the function Sfun[m,γ]  where K=Kg and Δx=deltax

Sfun[m_,gamma_]:=(Kg*Gamma[2-gamma])/(deltax)^2*(tt[m ]-tt[m-1])^gamma

Next the matrix system 

is built and then solved inside the functions  U1stepMatrixt  and U2stepsMatrixt.
The matrix A is built by means of the instruction SparseArray[...], while the vector G is generated by means of the instruction Table[Which[j==0,.....,{j,0,Np}] . Here uL[tt[m]]=u((Subscript[x, 0]),(Subscript[t, m]))  is the boundary condition to the left (at x=(Subscript[x, 0])) and uR[tt[m]]=u((Subscript[x, L]),(Subscript[t, m]))  is the boundary condition to the right (at x=(Subscript[x, R])). The solution {(Subscript[U, 0]^((m))),(Subscript[U, 1]^((m))),(Subscript[U, 2]^((m))),... (Subscript[U, Np]^((m)))}, is stored in the vector uVectortm

U1stepMatrixt[m]   provides the solution of 

i.e., it provides the solution {(Subscript[U, 0]^((m))),(Subscript[U, 1]^((m))),(Subscript[U, 2]^((m))),... (Subscript[U, Np]^((m)))}  at time (Subscript[t, m]) from the solutions {{(Subscript[U, 0]^((m-1))),(Subscript[U, 1]^((m-1))),(Subscript[U, 2]^((m-1))),... (Subscript[U, Np]^((m-1)))},{(Subscript[U, 0]^((m-2))),(Subscript[U, 1]^((m-2))),(Subscript[U, 2]^((m-2))),... (Subscript[U, Np]^((m-2)))} ,⋯} ,  employing a single timestep of size  (Subscript[t, m])-(Subscript[t, m-1])=ht

U1stepMatrixt[m_]:=(tt[0]=0;
tt[m]= tt[m-1]+ht;
uVectortm=
LinearSolve[
SparseArray[{Band[{1,2}]->Table[Which[j==0,0.,True,-Sfun[m,gammak]],{j,0,Np-1}],Band[{1,1}]->Table[Which[j==0,1.,j==Np,1.,True,1.+2.*Sfun[m,gammak]],{j,0,Np}],Band[{2,1}]->Table[Which[j==Np-1,0.,True,-Sfun[m,gammak]],{j,0,Np-1}]},Np+1],Table[Which[j==0,uL[tt[m]],j==Np,uR[tt[m]],True,GG[j,m,gammak]],{j,0,Np}]])

U2stepsMatrixt[n]  provides the solution {(Subscript[U, 0]^((n))),(Subscript[U, 1]^((n))),(Subscript[U, 2]^((n))),... (Subscript[U, Np]^((n)))}  at time (Subscript[t, n]) from the solutions {{(Subscript[U, 0]^((n-1))),(Subscript[U, 1]^((n-1))),(Subscript[U, 2]^((n-1))),... (Subscript[U, Np]^((n-1)))},{(Subscript[U, 0]^((n-2))),(Subscript[U, 1]^((n-2))),(Subscript[U, 2]^((n-2))),... (Subscript[U, Np]^((n-2)))} ,⋯} ,  employing  two timesteps of size  ((Subscript[t, n])-(Subscript[t, n-1]))/2=ht/2=hthalf

U2stepsMatrixt[n_]:=(
 hthalf=ht/2.;
Do[tt[m]= tt[m-1]+hthalf;
uVectortm=
LinearSolve[
SparseArray[{Band[{1,2}]->Table[Which[j==0,0.,True,-Sfun[m,gammak]],{j,0,Np-1}],Band[{1,1}]->Table[Which[j==0,1.,j==Np,1.,True,1.+2.*Sfun[m,gammak]],{j,0,Np}],Band[{2,1}]->Table[Which[j==Np-1,0.,True,-Sfun[m,gammak]],{j,0,Np-1}]},Np+1],Table[Which[j==0,uL[tt[m]],j==Np,uR[tt[m]],True,GG[j,m,gammak]],{j,0,Np}]];
Umatrix=Append[Umatrix,uVectortm],
{m,n,n+1}
]; 
Umatrix=Drop[Umatrix,-2];
uVectortm)

uMsolutionFixedTimesteps[t,Δ]  provides the full  solution   {(Subscript[U, 0]^((m))),(Subscript[U, 1]^((m))),(Subscript[U, 2]^((m))),... (Subscript[U, Np]^((m)))}  with m=0,1,2,....n  for fixed timesteps of size Δ from t=0 to t=(Subscript[t, n])  where n=(Subscript[t, n])/Δ .  The full solution is stored in the matrix Umatrix so that its element (1+m,1+j), i.e., Umatrix[[1+m,1+j]],  gives  (Subscript[U, j]^((m))).

uMsolutionFixedTimesteps[t_,htFixed_]:=Module[{numpasos},
ht=htFixed;
numpasos=Round[t/htFixed];
Umatrix={};Umatrix=AppendTo[Umatrix,UmatrixIni];
Do[
	Umatrix=Append[Umatrix,U1stepMatrixt[m]],
{m,1,numpasos}];
Umatrix
]

uMsolutionTE[n,τ]  provides the full  solution   {(Subscript[U, 0]^((m))),(Subscript[U, 1]^((m))),(Subscript[U, 2]^((m))),... (Subscript[U, Np]^((m)))}  with m=0,1,2,....n  for variable timesteps whose size is determined by the trial and error (T&E) method for a tolerance τ .  The full solution is stored in the matrix Umatrix so that its element (1+m,1+j), i.e., Umatrix[[1+m,1+j]],  gives  (Subscript[U, j]^((m))).

Some arbitrary parameters one has to chose:
• The initial  timestep (Subscript[Δ, 0])=Delta0TE.  Here is fixed to the arbitrary (but hardly relevant) value of 0.01
• The maximum value of the timestep, DeltatcapTE, one is ready to accept.  Here we choose Deltatmax=(10^(6 )) (the actual value could be somewhat larger). [When the solution goes towards zero, the absolute value of this solution will be smaller than the tolerance for times long enough. In this case, the condition =udif<tolerance always holds and the search of  ever-increasing timesteps would never end. Besides, to cap the size of the timesteps guarantees a minimum sampling of the solution.]

Delta0TE=0.01;DeltatcapTE=10.^6;

uMsolutionTE[numpasos_,tolerance_]:=(
Umatrix={};Umatrix=AppendTo[Umatrix,UmatrixIni];
ht=Delta0TE;
Do[
usingle=U1stepMatrixt[m];
uhalf=U2stepsMatrixt[m];
udif=Max[Abs[usingle-uhalf]];
If[udif<tolerance,
	While[	
	ht<DeltatcapTE&&udif<tolerance,
	ht=2.*ht;
	usingle=U1stepMatrixt[m];
	uhalf=U2stepsMatrixt[m];
	udif=Max[Abs[usingle-uhalf]];
		         ];
        ht= ht/2.;
    ,
	While[
	udif>tolerance,
	ht=ht/2;
	usingle=U1stepMatrixt[m];
	uhalf=U2stepsMatrixt[m];
	udif=Max[Abs[usingle-uhalf]];
	
	         ]
   ];
 usingle=U1stepMatrixt[m];
uhalf=U2stepsMatrixt[m];
udif=Max[Abs[usingle-uhalf]];

Umatrix=Append[Umatrix,U1stepMatrixt[m]],
{m,1,numpasos}];
Umatrix
)

uMsolutionPredictive[n,τ]  provides the full  solution   {(Subscript[U, 0]^((m))),(Subscript[U, 1]^((m))),(Subscript[U, 2]^((m))),... (Subscript[U, Np]^((m)))}  with m=0,1,2,....n  for variable timesteps whose size is determined by the predictive method for a tolerance τ .  The full solution is stored in the matrix Umatrix so that its element (1+m,1+j), i.e., Umatrix[[1+m,1+j]],  gives  (Subscript[U, j]^((m))).

Some arbitrary parameters one has to chose:
• The initial  timestep (Subscript[Δ, 0])=Delta0Predictive. Here is fixed to the arbitrary (but hardly relevant) value of 0.01
• The maximum value of the timestep, DeltatcapTE, one is ready to accept.  Here we choose Deltatmax=(10^(6 )) (the actual value could be somewhat larger). [When the solution goes towards zero, the absolute value of this solution will be smaller than the tolerance for times long enough. In this case, the condition =udif<fac2*tolerance always holds and the search of  ever-increasing timesteps would never end. Besides, to cap the size of the timesteps guarantees a minimum sampling of the solution.]
• The under-relaxation parameter ω=omegaParameter. Here we use  ω=1/2. 
• The scaling exponent θ=thetaParameter.  Here we use  θ=3/2. 
• The fac1 and fac2 parameteres. The timestep is accepted if the difference =udif  is larger than fac1*tolerance and smaller than fac2*tolerance.  Here we use fac1=1/2 and fac2=2  [see Eq. (18)]

Delta0Predictive=0.01;DeltatcapPredictive=10.^6;omegaParameter=0.5;thetaParameter=3./2.; fac1=0.5; fac2=2.;

uMsolutionPredictive[numpasos_,tolerance_]:=(
Umatrix={};Umatrix=AppendTo[Umatrix,UmatrixIni];
ht=Delta0Predictive; 
Do[
usingle=U1stepMatrixt[m];
uhalf=U2stepsMatrixt[m];
udif=Max[Abs[usingle-uhalf]];
	While[	ht<DeltatcapPredictive&&(udif<fac1*tolerance||udif>fac2*tolerance),	ht=(1.-omegaParameter)*ht+omegaParameter*ht*(tolerance/udif)^(1./thetaParameter);  
	uhalf=U2stepsMatrixt[m];
	usingle=U1stepMatrixt[m];	
	udif=Max[Abs[usingle-uhalf]];
		         ];
 
Umatrix=Append[Umatrix,U1stepMatrixt[m]],
{m,1,numpasos}];
htVariableUltimo=ht;
Umatrix
)

Example 1

Here we provide the parameters of the problem and its exact solution

Parameter:  (Subscript[K, γ])=Kg=1

Kg=1.;

Exact solution:

uSol[x_,t_,gamma_]:=Sin[x]*MittagLefflerE[gamma,1,-t^gamma]

uSol[x_,t_,1]:=Sin[x]*(E^(-t))

Source term.  In this example is zero.

funSource[x_,t_,ga_]:=0.

Order of the fractional derivative: γ=gammak. Here we choose γ=1/4

gammak=0.25;

Spatial discretization:   Δx≡deltax=Pi/Np where, in this case, we choose  Np=2Nx=40

Nx=20;Np=2*Nx;deltax=Pi/(2.*Nx);

We set the initial condition:  (U^((0)))=UmatrixIni

UmatrixIni=Table[ Sin[j*deltax] ,{j,0,Np}]//N

Boundary conditions:  u(x=0,t)=uL[t]=0;   u(x=π,t)=uR[t]=0.

uL[t_]=0.;
uR[t_]=0.;

■ Numerical results with the T&E method for tolerance  τ=0.001

uMs is a matrix that contains the full solution {(Subscript[U, 0]^((m))),(Subscript[U, 1]^((m))),(Subscript[U, 2]^((m))),... (Subscript[U, Np]^((m)))}  for m=0,1,2,....,numSteps.

tolerance=0.001;
numSteps=20;
CPUtime=Timing[uMs=uMsolutionTE[numSteps,tolerance]];
Print["Example 1. T&E method. Case with τ=",tolerance," for γ=",gammak,", Np=",Np,",\!\(TraditionalForm\`\(\(\\ \)\*SubscriptBox[\(Δ\), \(0\)]\)\)=",Delta0TE,", (Subscript[Δ, max])=" ,DeltatcapTE ];
Print["CPU time=",CPUtime[[1]]," seconds"];
Print["The solution is saved in the matrix uMs"];

The solution for the m-th timestep is just the m+1 element of this matrix:  uMs[[m+1]]={(Subscript[U, 0]^((m))),(Subscript[U, 1]^((m))),(Subscript[U, 2]^((m))),... (Subscript[U, Np]^((m)))}

mStep=5;
Print["The numerical solution at time ",tt[mStep]," is:"];uMs[[mStep]]

The [[m+1,j+1] element of this matrix is the numerical solution for time (Subscript[t, m]) and position (Subscript[x, j]):  uMs[[m+1,j+1]]=(Subscript[U, j]^((m)))

mStep=5;jPos=Nx+1;
Print["The numerical solution at time t=",tt[mStep]," and position x= " ,jPos*deltax ," is:"];uMs[[mStep,jPos]]

A plot comparing the numerical and the exact solution at time (Subscript[t, mStep])

mStep=5;
Print["T&E method"];
Print["Numerical solution (points) and exact solution (line)  at time t=",tt[mStep]];
 FigNum=ListPlot[Table[{xj=j*deltax,uMs[[mStep+1,j+1]]},{j,0,Np}],PlotRange->All];FigExact=Plot[uSol[x,tt[mStep],gammak],{x,0,Pi},PlotRange->All]; Show[FigNum,FigExact,Frame->True,FrameLabel->{"x","u"}]

Here we give a table that compares the numerical and exact solution at the middle point, x=π/2, for the first numSteps timesteps.

headTE={"step","t","(Subscript[U, T & E])(π/2,t)","(Subscript[u, exact])(π/2,t)","error"};

tabTE=Table[{m,tt[m],uMs[[m+1,Nx+1]],uSol[Pi/2,tt[m],gammak],uMs[[m+1,Nx+1]]-uSol[Pi/2,tt[m],gammak]},{m,0,numSteps}]//N;Insert[tabTE,headTE,1]//TableForm

A plot comparing the numerical and the exact solution at  position x=π/2 for times  t  with   0≤t≤ (Subscript[t, numSteps])

Print["T&E method"];Print["Numerical solution (points) and exact solution (line) at position x=π/2 vs. time"]; FigNum=ListPlot[Table[{tabTE[[n,2]],tabTE[[n,3]]},{n,1,numSteps+1}],PlotRange->All];
FigExact=Plot[uSol[Pi/2,time,gammak],{time,0,1.1*tabTE[[numSteps+1,2]]},PlotRange->All];
Show[FigNum,FigExact,Frame->True,FrameLabel->{"time","u"},PlotRange->All]

■ Numerical results with the predictive method for tolerance  τ=0.001

uMs is a matrix that contains the full solution {(Subscript[U, 0]^((m))),(Subscript[U, 1]^((m))),(Subscript[U, 2]^((m))),... (Subscript[U, Np]^((m)))}  for m=0,1,2,....,numSteps.

tolerance=0.001;
numSteps=20;
CPUtime=Timing[uMs=uMsolutionPredictive[numSteps,tolerance]];
Print["Example 1. Predictive method. Case with τ=",tolerance," for γ=",gammak,", Np=",Np,",\!\(TraditionalForm\`\(\(\\ \)\*SubscriptBox[\(Δ\), \(0\)]\)\)=",Delta0Predictive,", (Subscript[Δ, max])=" ,DeltatcapPredictive  , ",ω=",  omegaParameter, ",\!\(TraditionalForm\`\(\(\\ \)\(θ\)\)\)=",  thetaParameter  ];
Print["CPU time=",CPUtime[[1]]," seconds"];
Print["The solution is saved in the matrix uMs"];

The solution for the m-th timestep is just the m+1 element of this matrix:  uMs[[m+1]]={(Subscript[U, 0]^((m))),(Subscript[U, 1]^((m))),(Subscript[U, 2]^((m))),... (Subscript[U, Np]^((m)))}

mStep=5;
Print["The numerical solution at time ",tt[mStep]," is:"];uMs[[mStep]]

The [[m+1,j+1] element of this matrix is the numerical solution for time (Subscript[t, m]) and position (Subscript[x, j]):  uMs[[m+1,j+1]]=(Subscript[U, j]^((m)))

mStep=5;jPos=Nx+1;
Print["The numerical solution at time t=",tt[mStep]," and position x= " ,jPos*deltax ," is:"];uMs[[mStep,jPos]]

A plot comparing the numerical and the exact solution at time (Subscript[t, mStep])

mStep=5;
Print["Predictive method"];
Print["Numerical solution (points) and exact solution (line)  at time t=",tt[mStep]];
 FigNum=ListPlot[Table[{xj=j*deltax,uMs[[mStep+1,j+1]]},{j,0,Np}],PlotRange->All];FigExact=Plot[uSol[x,tt[mStep],gammak],{x,0,Pi},PlotRange->All]; Show[FigNum,FigExact,Frame->True,FrameLabel->{"x","u"}]

Here we give a table that compares the numerical and exact solution at the middle point, x=π/2, for the first numSteps timesteps.

headPredictive={"step","t","(Subscript[U, Predic])(π/2,t)","(Subscript[u, exact])(π/2,t)","error"};

tabPred=Table[{m,tt[m],uMs[[m+1,Nx+1]],uSol[Pi/2,tt[m],gammak],uMs[[m+1,Nx+1]]-uSol[Pi/2,tt[m],gammak]},{m,0,numSteps}]//N;Insert[tabPred,headPredictive,1]//TableForm

A plot comparing the numerical and the exact solution at  position x=π/2 for times  t  with   0≤t≤ (Subscript[t, numSteps])

Print["Predictive method"];Print["Numerical solution (points) and exact solution (line) at position x=π/2 vs. time"]; FigNum=ListPlot[Table[{tabPred[[n,2]],tabPred[[n,3]]},{n,1,numSteps}],PlotRange->All];
FigExact=Plot[uSol[Pi/2,time,gammak],{time,0,1.1*tabPred[[numSteps+1,2]]},PlotRange->All];
Show[FigNum,FigExact,Frame->True,FrameLabel->{"time","u"},PlotRange->All]

■ Numerical results with the method with fixed timesteps for Δ=0.01=htFT

uMs is a matrix that contains the full solution {(Subscript[U, 0]^((m))),(Subscript[U, 1]^((m))),(Subscript[U, 2]^((m))),... (Subscript[U, Np]^((m)))}  for m=0,1,2,....,numSteps.

htFT=0.01;
numSteps=50;
CPUtime=Timing[uMs=uMsolutionFixedTimesteps[numSteps*htFT,htFT]];
Print["Example 1. Method with fixed timesteps of size ",htFT," for γ=",gammak,", Np=",Np ];
Print["CPU time=",CPUtime[[1]]," seconds"];
Print["The solution is saved in the matrix uMs"];

The solution for the m-th timestep is just the m+1 element of this matrix:  uMs[[m+1]]={(Subscript[U, 0]^((m))),(Subscript[U, 1]^((m))),(Subscript[U, 2]^((m))),... (Subscript[U, Np]^((m)))}

mStep=5;
Print["The numerical solution at time ",tt[mStep]," is:"];uMs[[mStep]]

The [[m+1,j+1] element of this matrix is the numerical solution for time (Subscript[t, m]) and position (Subscript[x, j]):  uMs[[m+1,j+1]]=(Subscript[U, j]^((m)))

mStep=5;jPos=Nx+1;
Print["The numerical solution at time t=",tt[mStep]," and position x= " ,jPos*deltax ," is:"];uMs[[mStep,jPos]]

A plot comparing the numerical and the exact solution at time (Subscript[t, mStep])

mStep=5;
Print["Method with fixed timesteps"];
Print["Numerical solution (points) and exact solution (line)  at time t=",tt[mStep]];
 FigNum=ListPlot[Table[{xj=j*deltax,uMs[[mStep+1,j+1]]},{j,0,Np}],PlotRange->All];FigExact=Plot[uSol[x,tt[mStep],gammak],{x,0,Pi},PlotRange->All]; Show[FigNum,FigExact,Frame->True,FrameLabel->{"x","u"}]

Here we give a table that compares the numerical and exact solution at the middle point, x=π/2, for the first numSteps timesteps.

headFixedTimesteps={"step","t","\!\(\*SubscriptBox[\(U\), \(Fixed\\ Δ\)]\)(π/2,t)","(Subscript[u, exact])(π/2,t)","error"};

tabFT=Table[{m,tt[m],uMs[[m+1,Nx+1]],uSol[Pi/2,tt[m],gammak],uMs[[m+1,Nx+1]]-uSol[Pi/2,tt[m],gammak]},{m,0,numSteps}]//N;Insert[tabFT,headFixedTimesteps,1]//TableForm

A plot comparing the numerical and the exact solution at  position x=π/2 for times  t  with   0≤t≤ (Subscript[t, numSteps])

Print["Method with fixed timesteps"];Print["Numerical solution (points) and exact solution (line) at position x=π/2 vs. time"]; FigNum=ListPlot[Table[{tabFT[[n,2]],tabFT[[n,3]]},{n,1,numSteps}],PlotRange->All];
FigExact=Plot[uSol[Pi/2,time,gammak],{time,0,1.1*tabFT[[numSteps+1,2]]},PlotRange->All];
Show[FigNum,FigExact,Frame->True,FrameLabel->{"time","u"},PlotRange->All]

Example 2

Here we provide the parameters of the problem and its exact solution

Parameters:  (Subscript[K, γ])=Kg=1, ν=nu=2*Pi;  a=aParameter=1/10

Kg=1.;nu=2.*Pi;aParameter=0.1;

Exact solution:

uSol[x_,t_,ga_]:=Sin[x]*(MittagLefflerE[ga,-t^ga]+aParameter*Sin[nu*t])

Source term:

funSource[x_,t_,ga_]:=aParameter*(nu^ga*Sin[nu*t+ga*Pi/2]+Sin[nu*t])* Sin[x]

Order of the fractional derivative: γ=gammak. Here we choose γ=1/4

gammak=0.25;

Spatial discretization:   Δx≡deltax=Pi/Np where, in this case, we choose  Np=2Nx=40

Nx=20;Np=2*Nx;deltax=Pi/(2.*Nx);

We set the initial condition:  (U^((0)))=UmatrixIni

UmatrixIni=Table[ Sin[j*deltax] ,{j,0,Np}]//N

Boundary conditions:  u(x=0,t)=uL[t]=0;   u(x=π,t)=uR[t]=0.

uL[t_]=0.;
uR[t_]=0.;

■ Numerical results with the T&E method for tolerance  τ=0.001

uMs is a matrix that contains the full solution {(Subscript[U, 0]^((m))),(Subscript[U, 1]^((m))),(Subscript[U, 2]^((m))),... (Subscript[U, Np]^((m)))}  for m=0,1,2,....,numSteps.

tolerance=0.001;
numSteps=40;
CPUtime=Timing[uMs=uMsolutionTE[numSteps,tolerance]];
Print["Example 2. T&E method. Case with τ=",tolerance," for γ=",gammak,", Np=",Np,",\!\(TraditionalForm\`\(\(\\ \)\*SubscriptBox[\(Δ\), \(0\)]\)\)=",Delta0TE,", (Subscript[Δ, max])=" ,DeltatcapTE ];
Print["CPU time=",CPUtime[[1]]," seconds"];
Print["The solution is saved in the matrix uMs"];

The solution for the m-th timestep is just the m+1 element of this matrix:  uMs[[m+1]]={(Subscript[U, 0]^((m))),(Subscript[U, 1]^((m))),(Subscript[U, 2]^((m))),... (Subscript[U, Np]^((m)))}

mStep=5;
Print["The numerical solution at time ",tt[mStep]," is:"];uMs[[mStep]]

The [[m+1,j+1] element of this matrix is the numerical solution for time (Subscript[t, m]) and position (Subscript[x, j]):  uMs[[m+1,j+1]]=(Subscript[U, j]^((m)))

mStep=5;jPos=Nx+1;
Print["The numerical solution at time t=",tt[mStep]," and position x= " ,jPos*deltax ," is:"];uMs[[mStep,jPos]]

A plot comparing the numerical and the exact solution at time (Subscript[t, mStep])

mStep=5;
Print["T&E method"];
Print["Numerical solution (points) and exact solution (line)  at time t=",tt[mStep]];
 FigNum=ListPlot[Table[{xj=j*deltax,uMs[[mStep+1,j+1]]},{j,0,Np}],PlotRange->All];FigExact=Plot[uSol[x,tt[mStep],gammak],{x,0,Pi},PlotRange->All]; Show[FigNum,FigExact,Frame->True,FrameLabel->{"x","u"}]

Here we give a table that compares the numerical and exact solution at the middle point, x=π/2, for the first numSteps timesteps.

headTE={"step","t","(Subscript[U, T & E])(π/2,t)","(Subscript[u, exact])(π/2,t)","error"};

tabTE=Table[{m,tt[m],uMs[[m+1,Nx+1]],uSol[Pi/2,tt[m],gammak],uMs[[m+1,Nx+1]]-uSol[Pi/2,tt[m],gammak]},{m,0,numSteps}]//N;Insert[tabTE,headTE,1]//TableForm

A plot comparing the numerical and the exact solution at  position x=π/2 for times  t  with   0≤t≤ (Subscript[t, numSteps])

Print["T&E method"];Print["Numerical solution (points) and exact solution (line) at position x=π/2 vs. time"]; FigNum=ListPlot[Table[{tabTE[[n,2]],tabTE[[n,3]]},{n,1,numSteps}],PlotRange->All];
FigExact=Plot[uSol[Pi/2,time,gammak],{time,0,1.1*tabTE[[numSteps+1,2]]},PlotRange->All];
Show[FigNum,FigExact,Frame->True,FrameLabel->{"time","u"},PlotRange->All]

■ Numerical results with the predictive method for tolerance  τ=0.001

uMs is a matrix that contains the full solution {(Subscript[U, 0]^((m))),(Subscript[U, 1]^((m))),(Subscript[U, 2]^((m))),... (Subscript[U, Np]^((m)))}  for m=0,1,2,....,numSteps.

tolerance=0.001;
numSteps=40;
CPUtime=Timing[uMs=uMsolutionPredictive[numSteps,tolerance]];
Print["Example 2. Predictive method. Case with τ=",tolerance," for γ=",gammak,", Np=",Np,",\!\(TraditionalForm\`\(\(\\ \)\*SubscriptBox[\(Δ\), \(0\)]\)\)=",Delta0Predictive,", (Subscript[Δ, max])=" ,DeltatcapPredictive ,  ",ω=",  omegaParameter, ",\!\(TraditionalForm\`\(\(\\ \)\(θ\)\)\)=",  thetaParameter  ];
Print["CPU time=",CPUtime[[1]]," seconds"];
Print["The solution is saved in the matrix uMs"];

The solution for the m-th timestep is just the m+1 element of this matrix:  uMs[[m+1]]={(Subscript[U, 0]^((m))),(Subscript[U, 1]^((m))),(Subscript[U, 2]^((m))),... (Subscript[U, Np]^((m)))}

mStep=5;
Print["The numerical solution at time ",tt[mStep]," is:"];uMs[[mStep]]

The [[m+1,j+1] element of this matrix is the numerical solution for time (Subscript[t, m]) and position (Subscript[x, j]):  uMs[[m+1,j+1]]=(Subscript[U, j]^((m)))

mStep=5;jPos=Nx+1;
Print["The numerical solution at time t=",tt[mStep]," and position x= " ,jPos*deltax ," is:"];uMs[[mStep,jPos]]

A plot comparing the numerical and the exact solution at time (Subscript[t, mStep])

mStep=5;
Print["Predictive method"];
Print["Numerical solution (points) and exact solution (line)  at time t=",tt[mStep]];
 FigNum=ListPlot[Table[{xj=j*deltax,uMs[[mStep+1,j+1]]},{j,0,Np}],PlotRange->All];FigExact=Plot[uSol[x,tt[mStep],gammak],{x,0,Pi},PlotRange->All]; Show[FigNum,FigExact,Frame->True,FrameLabel->{"x","u"}]

Here we give a table that compares the numerical and exact solution at the middle point, x=π/2, for the first numSteps timesteps.

headPredictive={"step","t","(Subscript[U, Predic])(π/2,t)","(Subscript[u, exact])(π/2,t)","error"};

tabPred=Table[{m,tt[m],uMs[[m+1,Nx+1]],uSol[Pi/2,tt[m],gammak],uMs[[m+1,Nx+1]]-uSol[Pi/2,tt[m],gammak]},{m,0,numSteps}]//N;Insert[tabPred,headPredictive,1]//TableForm

A plot comparing the numerical and the exact solution at  position x=π/2 for times  t  with   0≤t≤ (Subscript[t, numSteps])

Print["Predictive method"];Print["Numerical solution (points) and exact solution (line) at position x=π/2 vs. time"]; FigNum=ListPlot[Table[{tabPred[[n,2]],tabPred[[n,3]]},{n,1,numSteps}],PlotRange->All];
FigExact=Plot[uSol[Pi/2,time,gammak],{time,0,1.1*tabPred[[numSteps+1,2]]},PlotRange->All];
Show[FigNum,FigExact,Frame->True,FrameLabel->{"time","u"},PlotRange->All]

■ Numerical results with the method with fixed timesteps for Δ=0.01=htFT

uMs is a matrix that contains the full solution {(Subscript[U, 0]^((m))),(Subscript[U, 1]^((m))),(Subscript[U, 2]^((m))),... (Subscript[U, Np]^((m)))}  for m=0,1,2,....,numSteps.

htFT=0.01;
numSteps=50;
CPUtime=Timing[uMs=uMsolutionFixedTimesteps[numSteps*htFT,htFT]];
Print["Example 1. Method with fixed timesteps of size ",htFT," for γ=",gammak,", Np=",Np ];
Print["CPU time=",CPUtime[[1]]," seconds"];
Print["The solution is saved in the matrix uMs"];

The solution for the m-th timestep is just the m+1 element of this matrix:  uMs[[m+1]]={(Subscript[U, 0]^((m))),(Subscript[U, 1]^((m))),(Subscript[U, 2]^((m))),... (Subscript[U, Np]^((m)))}

mStep=5;
Print["The numerical solution at time ",tt[mStep]," is:"];uMs[[mStep]]

The [[m+1,j+1] element of this matrix is the numerical solution for time (Subscript[t, m]) and position (Subscript[x, j]):  uMs[[m+1,j+1]]=(Subscript[U, j]^((m)))

mStep=5;jPos=Nx+1;
Print["The numerical solution at time t=",tt[mStep]," and position x= " ,jPos*deltax ," is:"];uMs[[mStep,jPos]]

A plot comparing the numerical and the exact solution at time (Subscript[t, mStep])

mStep=5;
Print["Method with fixed timesteps"];
Print["Numerical solution (points) and exact solution (line)  at time t=",tt[mStep]];
 FigNum=ListPlot[Table[{xj=j*deltax,uMs[[mStep+1,j+1]]},{j,0,Np}],PlotRange->All];FigExact=Plot[uSol[x,tt[mStep],gammak],{x,0,Pi},PlotRange->All]; Show[FigNum,FigExact,Frame->True,FrameLabel->{"x","u"}]

Here we give a table that compares the numerical and exact solution at the middle point, x=π/2, for the first numSteps timesteps.

headFixedTimesteps={"step","t","\!\(\*SubscriptBox[\(U\), \(Fixed\\ Δ\)]\)(π/2,t)","(Subscript[u, exact])(π/2,t)","error"};

tabFT=Table[{m,tt[m],uMs[[m+1,Nx+1]],uSol[Pi/2,tt[m],gammak],uMs[[m+1,Nx+1]]-uSol[Pi/2,tt[m],gammak]},{m,0,numSteps}]//N;Insert[tabFT,headFixedTimesteps,1]//TableForm

A plot comparing the numerical and the exact solution at  position x=π/2 for times  t  with   0≤t≤ (Subscript[t, numSteps])

Print["Method with fixed timesteps"];Print["Numerical solution (points) and exact solution (line) at position x=π/2 vs. time"]; FigNum=ListPlot[Table[{tabFT[[n,2]],tabFT[[n,3]]},{n,1,numSteps}],PlotRange->All];
FigExact=Plot[uSol[Pi/2,time,gammak],{time,0,1.1*tabFT[[numSteps+1,2]]},PlotRange->All];
Show[FigNum,FigExact,Frame->True,FrameLabel->{"time","u"},PlotRange->All]

Example 3

Here we provide the parameters of the problem and its exact solution

Parameters:  (Subscript[K, γ])=Kg=1, (Subscript[u, 0])=u0val=1,  L=LL=4,  M=Mvalue=8

Kg=1.;u0value=1.;LL=4.;Mvalue=8;

Exact solution [Eq. (26)]
 
 where z=x/(((K (t^γ)))^(1/2))  and    (Subscript[z, c]) = 2L / (((K (t^γ)))^(1/2))   and M→∞

Fox H functions  (Subscript[H, 10]^11)[z|({{1,γ/2}, {0,1}})] =H1011[z,γ]  for γ=1/4, 1/3, 1/2, 2/3, 3/4 and 1:

H1011[z_,1/4]:=-(1/(Gamma[(7/8)]))z HypergeometricPFQ[{(1/8)},{(1/4),(3/8),(1/2),(5/8),(3/4),(7/8),(9/8)},-((z^8)/16777216)]+(1/5040)(5040+(1/(Gamma[(3/4)]))2520 (z^2) HypergeometricPFQ[{(1/4)},{(3/8),(1/2),(5/8),(3/4),(7/8),(9/8),(5/4)},-((z^8)/16777216)]-(1/(Gamma[(5/8)]))840 (z^3) HypergeometricPFQ[{(3/8)},{(1/2),(5/8),(3/4),(7/8),(9/8),(5/4),(11/8)},-((z^8)/16777216)]+(1/(Sqrt[π]))210 (z^4) HypergeometricPFQ[{(1/2)},{(5/8),(3/4),(7/8),(9/8),(5/4),(11/8),(3/2)},-((z^8)/16777216)]-(1/(Gamma[(3/8)]))42 (z^5) HypergeometricPFQ[{(5/8)},{(3/4),(7/8),(9/8),(5/4),(11/8),(3/2),(13/8)},-((z^8)/16777216)]+(1/(Gamma[(1/4)]))7 (z^6) HypergeometricPFQ[{(3/4)},{(7/8),(9/8),(5/4),(11/8),(3/2),(13/8),(7/4)},-((z^8)/16777216)]-(1/(Gamma[(1/8)]))(z^7) HypergeometricPFQ[{(7/8)},{(9/8),(5/4),(11/8),(3/2),(13/8),(7/4),(15/8)},-((z^8)/16777216)])

H1011[z_,1/3]:=1+(1/120) z (-(1/(Gamma[(5/6)]))120 HypergeometricPFQ[{(1/6)},{(1/3),(1/2),(2/3),(5/6),(7/6)},-((z^6)/46656)]+z ((1/(Gamma[(2/3)]))60 HypergeometricPFQ[{(1/3)},{(1/2),(2/3),(5/6),(7/6),(4/3)},-((z^6)/46656)]+z (-(1/(Sqrt[π]))20 HypergeometricPFQ[{(1/2)},{(2/3),(5/6),(7/6),(4/3),(3/2)},-((z^6)/46656)]+z ((1/(Gamma[(1/3)]))5 HypergeometricPFQ[{(2/3)},{(5/6),(7/6),(4/3),(3/2),(5/3)},-((z^6)/46656)]-(1/(Gamma[(1/6)]))z HypergeometricPFQ[{(5/6)},{(7/6),(4/3),(3/2),(5/3),(11/6)},-((z^6)/46656)]))))

H1011[z_,1/2]:=1-(z HypergeometricPFQ[{(1/4)},{(1/2),(3/4),(5/4)},-((z^4)/256)])/Gamma[(3/4)]+(1/6) (z^2) ((1/(Sqrt[π]))3 HypergeometricPFQ[{(1/2)},{(3/4),(5/4),(3/2)},-((z^4)/256)]-(z HypergeometricPFQ[{(3/4)},{(5/4),(3/2),(7/4)},-((z^4)/256)])/Gamma[(1/4)])

H1011[z_,2/3]:=1-(z HypergeometricPFQ[{(1/3)},{(2/3),(4/3)},((z^3)/27)])/Gamma[(2/3)]+((z^2) HypergeometricPFQ[{(2/3)},{(4/3),(5/3)},((z^3)/27)])/(2 Gamma[(1/3)])

H1011[z_,3/4]:=-(1/(Gamma[(5/8)]))z HypergeometricPFQ[{(1/8),(11/24),(19/24)},{(1/4),(3/8),(1/2),(5/8),(3/4),(7/8),(9/8)},-((27 (z^8))/16777216)]+(1/5040)(5040+(1/(Gamma[(1/4)]))2520 (z^2) HypergeometricPFQ[{(1/4),(7/12),(11/12)},{(3/8),(1/2),(5/8),(3/4),(7/8),(9/8),(5/4)},-((27 (z^8))/16777216)]-(1/(Gamma[-(1/8)]))840 (z^3) HypergeometricPFQ[{(3/8),(17/24),(25/24)},{(1/2),(5/8),(3/4),(7/8),(9/8),(5/4),(11/8)},-((27 (z^8))/16777216)]-(1/(Sqrt[π]))105 (z^4) HypergeometricPFQ[{(1/2),(5/6),(7/6)},{(5/8),(3/4),(7/8),(9/8),(5/4),(11/8),(3/2)},-((27 (z^8))/16777216)]-(1/(Gamma[-(7/8)]))42 (z^5) HypergeometricPFQ[{(5/8),(23/24),(31/24)},{(3/4),(7/8),(9/8),(5/4),(11/8),(3/2),(13/8)},-((27 (z^8))/16777216)]+(1/(Gamma[-(5/4)]))7 (z^6) HypergeometricPFQ[{(3/4),(13/12),(17/12)},{(7/8),(9/8),(5/4),(11/8),(3/2),(13/8),(7/4)},-((27 (z^8))/16777216)]-(1/(Gamma[-(13/8)]))(z^7) HypergeometricPFQ[{(7/8),(29/24),(37/24)},{(9/8),(5/4),(11/8),(3/2),(13/8),(7/4),(15/8)},-((27 (z^8))/16777216)])

H1011[z_,1]:=Erfc[z/2]

Auxiliary functions for evaluating the exact solution

HEsca[z_,ga_]:=Which[z>50,0,z≤50,H1011[SetPrecision[z,60],ga]]    (*for z larger than 50 we approximate this function by zero*)

uEscalon[x_,t_,ga_,LL_,nterm_]:=(
z=x/(Kg*t)^(ga/2);
zc=(2*LL)/(Kg*t)^(ga/2);
HEsca[z,ga] -Sum[HEsca[n*zc-z,ga]-HEsca[n*zc+z,ga],{n,1,nterm}]
)  (*for very small times, so that z>50, this function is approximated by zero*)

Exact solution (an approximation retaining Mvalue terms in the series expansion)

uSol[x_,0,ga_]:=0/;x>0

uSol[x_,t_,ga_]:=u0value*uEscalon[x,t,ga,LL,Mvalue]

Source term.  In this example is zero.

funSource[x_,t_,ga_]:=0.

Order of the fractional derivative: γ=gammak. Here we choose γ=1/4

gammak=0.25;
gammaVal=Rationalize[gammak];

Spatial discretization:   Δx≡deltax=Pi/Np where, in this case, we choose  Np=2Nx=40

Nx=20;Np=2*Nx;deltax=LL/(2.*Nx);

We set the initial condition:  (U^((0)))=UmatrixIni

UmatrixIni=Table[ Which[j==0,u0value,j>0,0.] ,{j,0,Np}]

Boundary conditions:  u(x=0,t)=uL[t]=u0val=1;   u(x=π,t)=uR[t]=0.

uL[t_]=u0value;
uR[t_]=0.;

■ Numerical results with the T&E method for tolerance  τ=0.001

uMs is a matrix that contains the full solution {(Subscript[U, 0]^((m))),(Subscript[U, 1]^((m))),(Subscript[U, 2]^((m))),... (Subscript[U, Np]^((m)))}  for m=0,1,2,....,numSteps.

tolerance=0.001;
numSteps=50;
CPUtime=Timing[uMs=uMsolutionTE[numSteps,tolerance]];
Print["Example 3. T&E method. Case with Kg=",Kg,", (Subscript[u, 0])=",u0value,", L=",LL, ", M=",Mvalue];
Print["Numerical T&E solution for tolerance τ=",tolerance," γ=",gammak,", Np=",Np,",\!\(TraditionalForm\`\(\(\\ \)\*SubscriptBox[\(Δ\), \(0\)]\)\)=",Delta0TE,", (Subscript[Δ, max])=" ,DeltatcapTE ];
Print["CPU time=",CPUtime[[1]]," seconds"];
Print["The solution is saved in the matrix uMs"];

The solution for the m-th timestep is just the m+1 element of this matrix:  uMs[[m+1]]={(Subscript[U, 0]^((m))),(Subscript[U, 1]^((m))),(Subscript[U, 2]^((m))),... (Subscript[U, Np]^((m)))}

mStep=10;
Print["The numerical solution at time ",tt[mStep]," is:"];uMs[[mStep]]

The [[m+1,j+1] element of this matrix is the numerical solution for time (Subscript[t, m]) and position (Subscript[x, j]):  uMs[[m+1,j+1]]=(Subscript[U, j]^((m)))

mStep=10;jPos=Nx+1;
Print["The numerical solution at time t=",tt[mStep]," and position x= " ,jPos*deltax ," is:"];uMs[[mStep,jPos]]

A plot comparing the numerical and the exact solution at time (Subscript[t, mStep])

mStep=20;
Print["T&E method"];
Print["Numerical solution (points) and exact solution (line)  at time t=",tt[mStep]];
 FigNum=ListPlot[Table[{xj=j*deltax,uMs[[mStep+1,j+1]]},{j,0,Np}],PlotRange->All];FigExact=Plot[uSol[x,tt[mStep],gammaVal],{x,0,LL},PlotRange->All]; Show[FigNum,FigExact,Frame->True,FrameLabel->{"x","u"},PlotRange->All]

Here we give a table that compares the numerical and exact solution at the point xpos=1, for the first numSteps timesteps.

headTE={"step","t","(Subscript[U, T & E])(1,t)","(Subscript[u, exact])(1,t)","error"};

xPos=1.;
jPos=Round[xPos/deltax];tabTE=Table[{m,tt[m],uMs[[m+1,jPos+1]],uSol[xPos,tt[m],gammaVal],uMs[[m+1,jPos+1]]-uSol[xPos,tt[m],gammaVal]},{m,0,numSteps}]//N;Insert[tabTE,headTE,1]//TableForm

A plot comparing the numerical and the exact solution at  position x=xpos for times  t  with   0≤t≤ (Subscript[t, numSteps])

Print["T&E method"];Print["Numerical solution (points) and exact solution (line) at position x=",xPos," vs. time"]; FigNum=ListPlot[Table[{tabTE[[n,2]],tabTE[[n,3]]},{n,1,numSteps}],PlotRange->All];
FigExact=Plot[uSol[xPos,time,gammaVal],{time,0,1.1*tabTE[[numSteps+1,2]]},PlotRange->All];
Show[FigNum,FigExact,Frame->True,FrameLabel->{"time","u"},PlotRange->All]

■ Numerical results with the predictive method for tolerance  τ=0.001

uMs is a matrix that contains the full solution {(Subscript[U, 0]^((m))),(Subscript[U, 1]^((m))),(Subscript[U, 2]^((m))),... (Subscript[U, Np]^((m)))}  for m=0,1,2,....,numSteps.

tolerance=0.001;
numSteps=50;
CPUtime=Timing[uMs=uMsolutionPredictive[numSteps,tolerance]];
Print["Example 3. Predictive method. Case with Kg=",Kg,", (Subscript[u, 0])=",u0value,", L=",LL, ", M=",Mvalue];
Print["Numerical Predictive solution for tolerance τ=",tolerance," γ=",gammak,", Np=",Np,",\!\(TraditionalForm\`\(\(\\ \)\*SubscriptBox[\(Δ\), \(0\)]\)\)=",Delta0Predictive,", (Subscript[Δ, max])=" ,DeltatcapPredictive   ",ω=",  omegaParameter, ",\!\(TraditionalForm\`\(\(\\ \)\(θ\)\)\)=",  thetaParameter   ];
Print["CPU time=",CPUtime[[1]]," seconds"];
Print["The solution is saved in the matrix uMs"];

The solution for the m-th timestep is just the m+1 element of this matrix:  uMs[[m+1]]={(Subscript[U, 0]^((m))),(Subscript[U, 1]^((m))),(Subscript[U, 2]^((m))),... (Subscript[U, Np]^((m)))}

mStep=5;
Print["The numerical solution at time ",tt[mStep]," is:"];uMs[[mStep]]

The [[m+1,j+1] element of this matrix is the numerical solution for time (Subscript[t, m]) and position (Subscript[x, j]):  uMs[[m+1,j+1]]=(Subscript[U, j]^((m)))

mStep=5;jPos=Nx+1;
Print["The numerical solution at time t=",tt[mStep]," and position x= " ,jPos*deltax ," is:"];uMs[[mStep,jPos]]

A plot comparing the numerical and the exact solution at time (Subscript[t, mStep])

mStep=30;
Print["Predictive method"];
Print["Numerical solution (points) and exact solution (line)  at time t=",tt[mStep]];
 FigNum=ListPlot[Table[{xj=j*deltax,uMs[[mStep+1,j+1]]},{j,0,Np}],PlotRange->All];FigExact=Plot[uSol[x,tt[mStep],gammaVal],{x,0,LL},PlotRange->All]; Show[FigNum,FigExact,Frame->True,FrameLabel->{"x","u"},PlotRange->All]

Here we give a table that compares the numerical and exact solution at the point xpos=1, for the first numSteps timesteps.

headPredictive={"step","t","(Subscript[U, Predic])(1,t)","(Subscript[u, exact])(1,t)","error"};

xPos=1.;
jPos=Round[xPos/deltax];tabPred=Table[{m,tt[m],uMs[[m+1,jPos+1]],uSol[xPos,tt[m],gammaVal],uMs[[m+1,jPos+1]]-uSol[xPos,tt[m],gammaVal]},{m,0,numSteps}]//N;Insert[tabPred,headPredictive,1]//TableForm

A plot comparing the numerical and the exact solution at  position x=xpos for times  t  with   0≤t≤ (Subscript[t, numSteps])

Print["Predictive method"];Print["Numerical solution (points) and exact solution (line) at position x=",xPos," vs. time"]; FigNum=ListPlot[Table[{tabPred[[n,2]],tabPred[[n,3]]},{n,1,numSteps}],PlotRange->All];
FigExact=Plot[uSol[xPos,time,gammaVal],{time,0,1.1*tabPred[[numSteps+1,2]]},PlotRange->All];
Show[FigNum,FigExact,Frame->True,FrameLabel->{"time","u"},PlotRange->All]

■ Numerical results with the method with fixed timesteps for Δ=0.01=htFT

uMs is a matrix that contains the full solution {(Subscript[U, 0]^((m))),(Subscript[U, 1]^((m))),(Subscript[U, 2]^((m))),... (Subscript[U, Np]^((m)))}  for m=0,1,2,....,numSteps.

htFT=0.01;
numSteps=50;
CPUtime=Timing[uMs=uMsolutionFixedTimesteps[numSteps*htFT,htFT]];
Print["Example 3. Method with fixed timesteps. Case with Kg=",Kg,", (Subscript[u, 0])=",u0value,", L=",LL];
Print["Numerical solution with fixed timesteps of size ",htFT," for γ=",gammak,", Np=",Np  ];
Print["CPU time=",CPUtime[[1]]," seconds"];
Print["The solution is saved in the matrix uMs"];

The solution for the m-th timestep is just the m+1 element of this matrix:  uMs[[m+1]]={(Subscript[U, 0]^((m))),(Subscript[U, 1]^((m))),(Subscript[U, 2]^((m))),... (Subscript[U, Np]^((m)))}

mStep=5;
Print["The numerical solution at time ",tt[mStep]," is:"];uMs[[mStep]]

The [[m+1,j+1] element of this matrix is the numerical solution for time (Subscript[t, m]) and position (Subscript[x, j]):  uMs[[m+1,j+1]]=(Subscript[U, j]^((m)))

mStep=5;jPos=Nx+1;
Print["The numerical solution at time t=",tt[mStep]," and position x= " ,jPos*deltax ," is:"];uMs[[mStep,jPos]]

A plot comparing the numerical and the exact solution at time (Subscript[t, mStep])

mStep=5;
Print["Method with fixed timesteps"];
Print["Numerical solution (points) and exact solution (line)  at time t=",tt[mStep]];
 FigNum=ListPlot[Table[{xj=j*deltax,uMs[[mStep+1,j+1]]},{j,0,Np}],PlotRange->All];FigExact=Plot[uSol[x,tt[mStep],gammaVal],{x,0,LL},PlotRange->All]; Show[FigNum,FigExact,Frame->True,FrameLabel->{"x","u"}]

Here we give a table that compares the numerical and exact solution at the middle point, x=xPos, for the first numSteps timesteps.

headFixedTimesteps={"step","t","\!\(\*SubscriptBox[\(U\), \(Fixed\\ Δ\)]\)(xPos,t)","(Subscript[u, exact])(xPos,t)","error"};

xPos=1.;
jPos=Round[xPos/deltax];
Print["xPos=",xPos];tabFT=Table[{m,tt[m],uMs[[m+1,jPos+1]],uSol[xPos,tt[m],gammaVal],uMs[[m+1,jPos+1]]-uSol[xPos,tt[m],gammaVal]},{m,0,numSteps}]//N;Insert[tabFT,headFixedTimesteps,1]//TableForm

A plot comparing the numerical and the exact solution at  position x=xPos for times  t  with   0≤t≤ (Subscript[t, numSteps])

Print["Method with fixed timesteps"];Print["Numerical solution (points) and exact solution (line) at position x=",xPos," vs. time"]; FigNum=ListPlot[Table[{tabFT[[n,2]],tabFT[[n,3]]},{n,1,numSteps}],PlotRange->All];
FigExact=Plot[uSol[xPos,time,gammaVal],{time,0,1.1*tabFT[[numSteps+1,2]]},PlotRange->All];
Show[FigNum,FigExact,Frame->True,FrameLabel->{"time","u"},PlotRange->All]